# Paper: Is it still the economy? Economic voting in polarized politics
# Individual-level simulation
# Replication: Figure S12 (Supplementary Material)
remove(list = ls())

library(conflicted)
library(tidyverse)
library(dplyr)
conflict_prefer("filter","dplyr")
library(ggplot2)
library(here) # When we have an R-project, it reads files without setting an directory
conflict_prefer("here", "here")
library(nnet)
library(ggeffects)
library(MASS)
conflict_prefer("select","dplyr")

load(here("Final_Paper", "R&R", "IndLevel.RData"))

table(model_data$vote_choice)
table(model_data$vote_inc_choice)

model_data <- model_data %>% 
  mutate(vote_inc_choice4 = ifelse(vote_choice == "democrat" & year == 1996, 
                                   "Inc", NA),
         vote_inc_choice4 = ifelse(vote_choice == "democrat" & year == 2000, 
                                   "Inc", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "democrat" & year == 2012, 
                                   "Inc", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "democrat" & year == 2016, 
                                   "Inc", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "republican" & year == 1992, 
                                   "Inc", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "republican" & year == 2004, 
                                   "Inc", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "republican" & year == 2008, 
                                   "Inc", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "republican" & year == 2020, 
                                   "Inc", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "republican" & year == 1996, 
                                   "Ops", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "republican" & year == 2000, 
                                   "Ops", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "republican" & year == 2012, 
                                   "Ops", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "republican" & year == 2016, 
                                   "Ops", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "democrat" & year == 1992, 
                                   "Ops", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "democrat" & year == 2004, 
                                   "Ops", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "democrat" & year == 2008, 
                                   "Ops", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "democrat" & year == 2020, 
                                   "Ops", vote_inc_choice4),
         vote_inc_choice4 = ifelse(vote_choice == "abs", 
                                   "Abs", vote_inc_choice4))

model_data <- model_data %>% 
  mutate(vote_inc_choice4 = ifelse(vote_choice == "other", 
                                   "other", vote_inc_choice4)) %>% 
  mutate(vote_inc_choice4 = factor(vote_inc_choice4, levels = c("Inc",
                                                                "Ops",
                                                                "Abs",
                                                                "other")))


inpart <- model_data %>% 
  filter(partisan == "Inpart")

summary(ml_inp <- multinom(vote_inc_choice4 ~ econ_eval*percep_polar + 
                             female + black + hispanic + asian + other_race + 
                             log_age + college + income + state_fips + 
                             election, data = inpart))

ml_inp2 <- ml_inp %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_inp2)
ml_inp2 <- data.frame(ml_inp2)

ml_inp2 <- ml_inp2 %>% 
  mutate(partisan = "Inpartisan")

outpart <- model_data %>% 
  filter(partisan == "Outpart")

summary(ml_out <- multinom(vote_inc_choice4 ~ econ_eval*percep_polar + 
                             female + black + hispanic + asian + other_race + 
                             log_age + college + income + state_fips + 
                             election, data = outpart))

ml_out2 <- ml_out %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_out2)
ml_out2 <- data.frame(ml_out2)

ml_out2 <- ml_out2 %>% 
  mutate(partisan = "Outpartisan")


indpart <- model_data %>% 
  filter(partisan == "Ind")

summary(ml_ind <- multinom(vote_inc_choice4 ~ econ_eval*percep_polar + 
                             female + black + hispanic + asian + other_race + 
                             log_age + college + income + state_fips + 
                             election, data = indpart))

ml_ind2 <- ml_ind %>% 
  ggeffect(terms = c("econ_eval", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_ind2)
ml_ind2 <- data.frame(ml_ind2)

ml_ind2 <- ml_ind2 %>% 
  mutate(partisan = "Independent")


ml_data <- rbind(ml_inp2, ml_out2, ml_ind2)
glimpse(ml_data)

g_account <- ml_data %>% 
  mutate(partisan = factor(partisan,
                           levels = c("Inpartisan", 
                                      "Independent",
                                      "Outpartisan"),
                           labels = c("President's Party", 
                                      "Independent",
                                      "Opposing Party"))) %>%  
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(x == 1) %>% 
  mutate(x = as.factor(x)) %>%
  filter(response.level == "Inc") %>% 
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(shape = group)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = partisan, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85, position = position_dodge(width = .5),
                  color = "darkslategray") + 
  theme_light() + 
  labs(x = NULL, y = "Probability of Voting for the Incumbent Party",
       title = "(a) Voting for the Incumbent Party") +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", "0.5", "1"),
                     values = c(15, 16, 17)) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 11),
        legend.position = "none",
        legend.direction = "horizontal",
        plot.title = element_text(size = 15))
g_account

g_opp <- ml_data %>% 
  mutate(partisan = factor(partisan,
                           levels = c("Inpartisan", 
                                      "Independent",
                                      "Outpartisan"),
                           labels = c("President's Party", 
                                      "Independent",
                                      "Opposing Party"))) %>%  
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(x == 1) %>% 
  mutate(x = as.factor(x)) %>%
  filter(response.level == "Ops") %>% 
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(shape = group)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = partisan, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85, position = position_dodge(width = .5),
                  color = "darkslategray") + 
  theme_light() + 
  labs(x = NULL, y = "Probability of Voting for the Opposition",
       title = "(b) Voting for the Opposition") +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", "0.5", "1"),
                     values = c(15, 16, 17)) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 11),
        legend.position = "none",
        legend.direction = "horizontal",
        plot.title = element_text(size = 15))
g_opp


g_abs <- ml_data %>% 
  mutate(partisan = factor(partisan,
                           levels = c("Inpartisan", 
                                      "Independent",
                                      "Outpartisan"),
                           labels = c("President's Party", 
                                      "Independent",
                                      "Opposing Party"))) %>%  
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(x == 1) %>% 
  mutate(x = as.factor(x)) %>%
  filter(response.level == "Abs") %>% 
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(shape = group)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = partisan, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85, position = position_dodge(width = .5),
                  color = "darkslategray") + 
  theme_light() + 
  labs(x = "Partisan Attachment", y = "Probability of Abstention",
       title = "(c) Abstention") +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", "0.5", "1"),
                     values = c(15, 16, 17)) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 12),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 15))
g_abs

glimpse(ml_data)
table(ml_data$response.level)
g_oth <- ml_data %>% 
  mutate(partisan = factor(partisan,
                           levels = c("Inpartisan", 
                                      "Independent",
                                      "Outpartisan"),
                           labels = c("President's Party", 
                                      "Independent",
                                      "Opposing Party"))) %>%  
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(x == 1) %>% 
  mutate(x = as.factor(x)) %>%
  filter(response.level == "other") %>% 
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(shape = group)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = partisan, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85, position = position_dodge(width = .5),
                  color = "darkslategray") + 
  theme_light() + 
  labs(x = "Partisan Attachment", y = "Probability of Voting for a 3rd Party",
       title = "(d) Voting for a Thirdy Party") +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", "0.5", "1"),
                     values = c(15, 16, 17)) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 12),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 15))
g_oth

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

multiplot(g_account, g_abs, g_opp, g_oth, cols = 2)

